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ABSTRACT 

This paper presents a multiphase turbulence 
closure employing one transport equation, namely 
the turbulence kinetic energy equation. Hie 
proposed form of this equation Is different from 
the earlier formulations in some aspects Ihe 

power spectrum of the carrier fluid is divided 
into two regions, which Interact in different ways 
and at different rates with the suspended 
particles as a function of the particle-eddy size 
ratio and density ratio. The length scale is 
described algebraically. A mass/time averaging 
procedure for the momentum and kinetic energy 
equations is adopted. The resulting turbulence 
correlations are modeled under less restrictive 
assumptions comparative to the previous work. The 
closures for the momentum and kinetic energy- 
equations are given. Comparisons of the 

predictions with experimental results on 
liquid-solid jet and gas-solid pipe flow show 
satisfactory agreement. 

NOMENCLATURE 


a amplitude ratio 

b body force 

r r C C (T, C = constants 

S>2* S3’ >r Ss’ nr m 

particle diameter 

turbulent diffusion coefficient for solid 
phase 

energy spectrum 
flow variable 
interaction term 

flux vector for a variable 4 


c r V 

d 

D 

E(»c) 

J. 


k kinetic energy of turbulence 

t length scale 

N Stokes number 

S 

K(r.t) phase distribution function, Eq (1) 

p pressure 

r posi t ion vector 

Re Reynolds number based on the most 

energetic eddy size 

t.AT.At time and averaging time intervals 
corresponding to the turbulence 


production and transfer range, 
respectively 

u,v Eulerian and Lagrangian velocities, 

respectively 

V.Av volume and averaging volume, respectively 
xxx Cartesian coordinates 

i ’ j ’ n 

y distance from the wall 

a concentration of volume 

t dissipation rate of turbulence kinetic 


energy 

Kolmogorov length scale 
wave number 
dynamic viscosity 
densi ty 

turbulent Schmidt/Prandt 1 number 
shear stress 


Subscripts 

i.j.n denote Cartesian coordinates ( = 1.2.3) 

e eddy 

k turbulence kinetic energy 

K flow component K 

KP phase K in the production range 

KT phase K in the transfer range 

£ laminar 

L,S denote liquid and solid, respectively 

p production range 

t turbulent 

T transfer range 

Tot total 

2 average over area 

t dissipation rate of turbulence kinetic 

energy 

Superscripts 


<f„> 


mass average of f 

turbulence fluctuation of f 

fluctuation of f at low wave num er 

(production range) 

fluctuation of f at inter 

number (transfer range) 

produc t ion 

transfer 

volume/time and mass/time averages of f 

over K r , v 

Intrinsic space average of f over h 


( ) 
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1- INTRODIK^TION 

Multiphase flows are widely applied in 
peering processes from chemical, petroleum 
mining and other industries. Various theoretical 

er\^r r technStpjes for the investigation 
°L7 ‘ Tr Pliable. Some of them are a 

straightforward extension from the single phase 
flow models by introducing some ad hoc 
modifications. Other investigations originate 

^d m mode'ls Ka f r° lld f !° W (S °°- 1983 > 0r fluidized 
oea models (Wang et al . . 19SS). 

Increasing concern for the prediction of 

durw7h r ltIphase fl °« have been noticed 
during the last twenty years (Danon et al. (1974) 

( 98m 7, V^ d3U ( 19773 • Genchev “ d Karpuzo^ 

W980) , Melville and Bray (1979) Crowe ,nH 

^iarna (1978). Michael ides and Farmer’ (1983) and 

Su'L k."®”' T "° *^ uatlon turbulence 

models have been proposed for dilute particulate 

Crowd ^ Elg . hob f shl et al - (1982. 1983. 1984) and 
Crowder et al - (1984). Algebraic and one equation 
turbulence models have been suggested also for 

1986) in* hlT k fl ° WS (R0C ° *' al " ,983> 1985 ’ 
i Uv - t h h thC partlcle- P a rticle interactions 

f Cid7n r ° besIdes the Huid-fluid and 

« " teractl ° ns - * ost of these studies 
litnitnHn ° V *“ r ll er Investigations have some 
limitations. In the above mentioned studies the 
response of solids to the turbulent fluctuations 

restrict! "VY, f,U,d Is ° b ^l"od under 
fl97S «° nS ^t m U v l ° lh ° Se refered W Hinze 

iddlri P ' 46 ° ) ‘ " hlch 1,mit thelr use. In 
add it on to that. empirical constants and 

empirical functions are usually introduced in 
these models. 

The purpose of the present paper Is 

T° P r °P<>« a specific mass/time averaging 
pproach for multiphase turbulent flows. Even if 
the approach is developed for incompressible flow, 
its application for other multiphase flows is 
foreseen: Proa the liquid-solid interaction forces 
only the drag force la considered in this peper. 

' To Improve the one-equation turbulence 
LIT!' I L C28] by lncludIn « ^e modulation 

»rtl7£ 7 ^ y J P artlcl « «* a function of 
particle size and density. 

mr.d .1 U1 J T ° test the proposed model with other 
modeis and experimental data for various two-phase 

coefficients 001 ad ° ptlng ^ adjusting empirical 
2. MIXED AVERAGING APPROACH 

continuum transport equations for 
multiphase flows can be obtained by assuming a 
continuum medium with averaged field qualities by 
using either time, local volume, local nass or 
nofi- 7 \ al * vera «ing (see Buyevich (1971). Soo 
J,™ l\l Ver T r " d Delta S re (1968). Hetsroni 
The averaging for multiphase flow 
systems may be performed in various ways. Mass 

(l«l lng / teCh " 1<, " t ras a PP lled W Abou-Arab 
' or turbulent Incompressible and 
nonun!f S b ! C i flows. To express the spatial 
ormities and interactions between the flow 

^ C ° and Sh ° 0k ( 1985 ) teve developed a 
turbuflnr V ° , . U T /tlme ave raging technique for 
slwif k™ 0mP ° ntnt systems. in which the 
size of the averaging volume Av is related to the 


of < he Eulerlaj> description 

description and there^are"' HT ^ U *l an S ian 
mathematical schemes 7 u comprehensive 

transformed the Volume/ tinJ* 0 ^ formulatlon - the y 
time averaging rl lm ? av e raging into double 

component^ ? n the ^ n ° W functlon f ^ d “>' 

K in the ■l J tture. at a position r and 
time t ~ 


>.t) = _L f 1 

AT J t . 


t+AT /2 


< f K > dt 


-here 1 

the turbu? 8108 * ntarv al corresponding to 
turbulence production range (4T ») 

< f > . _L f t+4t/2 

K ~ At J f (r.T) K(r ,r) R(r-t) dr 

J t-At/2 ‘ 

TTT* iC avera « ln S of f over At^d 
tne flow component K. 

At = Eulerlan time scale for the most energetic 

eddies =j" R( T -t)dT (corresponding to the 
o 

but 1 ™ k J ength SCale) - Note thac At«AT. 
at r H largCr than P ar tlcle residence time 

1 If the considered flow component 
K(£.t) = / resides at point r and time t , 

1 0 otherwise “ (2b) 

- phase distribution function 

R(r-t) = v "(t)Vft— t)~ 

(2c) 

V" 2 (t) 

= auto-correlation coefficient of the 
velocity fluctuations of the most 
energetic eddies (Taylor s scale) 

aver aging over At corresponds to the 
Av V0 } Ume Av ' ^ dimension of 

At = ni^ rh b 7 Cu / bulence D,ixin 6 length. and 
By lnt^^tinT^ver . 1V n ( r velocity) [29] 
Interfaces iuh other n C0mp0nenl K its 

5 rs?i 

jj re a a traighforward manner in the 

S ™^.‘SS* - V f 7- » (■) - 

* *-*■>— *? 

and f £ : K 

'k ■ r K * • \ * f K • rj 


*K ~ f K ~ ( f K y 

= spatial nonuniformities within Av (or 
temporal nonuniformities within At) 

f K = < > - f* 
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— temporal nonuniformities of ( f^ ) within AT. 

Since the averaging domain (Av or At) has the 
dimensions of the mixing length or Eulerian time 
scale, respectively, the turbulence fluctuation f£ 

corresponds to the turbulence transfer range. The 
temporal nonuniformities f^ reflect the turbulence 

fluctuations in the production range. 

By averaging with formula (1) the point 
instantaneous conservation equations, one obtains 
the double time averaged equations. The 

formulation Is equivalent to the volume/time 
averaging. The momentum and kinetic energy 

equations are given in Appendix A. 

The local mass average of f over a flow 
component K is denoted f^. It is obtained by 

applying (1). in which the phase distribution 
function K(r.T) is weighted by the specific mass 

V 

In this paper we model the phase interaction 
by using spectral analysis and suggest a closure 
of the mass averaged equations for linear momentum 
and kinetic energy. The averaged equations are 
initially written with all the terms, and then 
simplified formulations for various flow 

conditions are suggested. 

3. ENERGY SPECTRUM AND SOLIDS - EDDY INTERACTION 

It is well accepted that turbulence is 
characterized by fluctuating motions defined by 
an energy spectrum (Tennekes and Lumely (1972)). 
Single time scale models, which are normally used 
for the prediction of turbulent flows. seems 
simplistic because different turbulent 

interactions are associated with different parts 
of the energy spectrum (Hanjalic’ et al . (1979)). 
A typical energy spectrum can be divided into 
three regions. The first region is the production 
region of large eddies and low wave number. The 
third region Is the dissipation region with small 
eddies and high wave number, in which the total 
kinetic energy produced at the lower wave number 

is dissipated. The intermediate range of wave 
numbers represents the Taylor’s transfer range. 
The total kinetic energy k of turbulence may be 
divided into production range (k p ) and transfer 

range (k^.) because there is negligible kinetic 

energy in the dissipation range: 



k 

~ k P + k T 

(6) 

where 




k p = 

1 

= 2 

~2 
u i • 

(7a) 

k T : 

1 

= 2 

"2 

u i 

(7b) 

V 

U i 

= fluctuating velocities in the 



production and transfer range, 
respectively . 



Figure 1. Schenatic showing the relative particle 
size to different eddy sizes in the 
energy spectrum. 


This partitioning of the energy spectrum was shown 
to be important for swirling flows (Chen (1986)). 
and hetrogeneous mixture flows such as two-phase 
jet (Al-Taweel and Landau (1977)). 

By using spectral analysis in conjunction 
with nass/time averaging some additional turbulent 
correlations will result in the mixture flow 
equations comparative to homogenous flows. These 
correlations can be classified into £ L Y g 

categories : 

i) Eddy-eddy interaction 
11) Eddy-mean flow interaction 

iii) Eddy-particle interaction 

iv) Particle-mean flow interaction 

v) Particle-particle interaction (for dense 
suspension flow). 

These correlations have to be modeled. Since the 
suspended particles may be of different sizes and 
different materials. their response to the carrier 
fluid fluctuations will vary as a function of the 
mean and fluctuating properties of the flow. The 
present work will consider a two-way — Interaction 
mechanism between solid particles and fluid 
vorticies in dilute suspensions. This interaction 
mechanism depends on the ratio between the 
particle size d g and the turbulent vortex (eddy) 

size 1 . These length scales are compared with 
the turbulence dissipation micro-scale (q) • 

To determine the particle-eddy — Interactio n 
the energy spectrum for multiphase flow system is 
divided into three typical zones (Figure 1)' 


"Large vortex zone” (#!)* where the 
turbulence energy is extracted from the mean 
flow by low frequency eddies. Here, the eddy 
length scale is larger than the particle 

size 
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*el > d S > ^ 


( 8 a) 


•'Medium vortex zone" (82). where the solid 
particles are about the same size with the 
vortex size , 1 . e . 


'e2 S d S 


> V 


(6b) 


•■Small vortex zone" (83). which would 
correspond to the Kolmogorov’s length scale. 


d S > J e3 


(8c) 


In zone 81 the solid particles generally follow 
the motion within a vortex, and have an energy 
dissipation effect. The particle response to the 
turbulent fluctuations (turbulence modulation) is 
fully determined (see Hinze. 1975). In the snail 
vortex zone 83 the solid particles can not 
significant ly a f fe ct the turbulence 

micros true ture . For the intermediate zone 82 a 
linear variation of the particle response is 
considered. This partitioning allows for the 
par ticles-nonunl form size eddies interaction to be 
efficiently modeled. 

The present closure formulation originates 
rrom the idea of subgrid scale modeling. If this 
idea is to be accepted, any flow quantity u. v. a. 

k * etc - ** separated into three parts 

according to (3). where and f” define the 

fluctuations in the production and transfer ranges 
of the energy spectrum, respectively. By starting 
from the particle equation of motion in its 
general form, the relation between the particle 
motion and different fluid eddies can be 
determined. and from here the fluctuation 

components f£ and f ” (see section 6). 

4. COMPOSED AVERAGED EQUATIONS 
4.1 Mean Flow Governing Equations 

The mass/time averaged momentum equation 
(Appendix A. Eq. A2) with f * = f • ♦ f" yields- 


a 


P K dt ( “K u Ki + VV + “K U K + a u Ki + “'“Ki ) 
Time rate change of the mean flow: connection 

4p K sr CVKi^j] 

J J 

Mean flam cau>ection 


+ p k [ V K i u Kj + Vki-Kj] 

Inertial effect 


+P K (ta K U Ki u Kj + VKi u KjJ + “KjC^Ki +a X u K ^ 
Collisional/Inertial Effects 

+ “Ki [ VKj + + ^ + S^ Kj ] 


+ U Ki °K U KJ 4 “Kj^i + °x u Ki u Kj + a K u Ki U K j 

+ °K U Ki U Ki + a K u Ki u Kj + ^ u Kj u Ki + Vk. u K^ 

= P K [ \ b K, 4 ^ + ^ + ^ 
force 

_ _3_ r - - — _ 

la, < p t< 4 4 “k p k + “k P K + 

Pressure effect 


. ^K <3o v aa,; da' 

[p K dx 4 4 Pk ST + P^ dT + ^ -...] 

1 i i 

4 S/Yk + VkiJi + ^^7, 

Frictional effect 

4 ^ T K2JI 4 ^iaji] + (9) 

Interaction 

Where K is a flow component, b^ Is the body force 
‘n the 1-th direction. and is the 

projection 


ln t * ie direction of the 

Interaction vector ( I ,) 

*~K -K 

related^To' t*h (9) , Contains correlations which are 
ranees of rh produ ction and transfer wave number 

mixed turbu ence spectrum. It contains 

aiso mixed correlations, 

4 2 ~ e k inetic Energy Eouatirm 


the mean ° W governing equation 

k^LtTc IZ aVCragln « form of the turbulent 
kinetic energy equation can be derived for a flow 

component K. The exact form of this equation for 

(eT ^AsT 16 T UrbU ' ent flow ls « lven 1" Appendix A 
AS) I, contains more than one hundred 
correlations which are related to the eddies in 

m^xed r r i° n transfer ranges as well as some 

are oredom' 13 * ' However ’ onl y some of these 
function of" 311 '. 8 giVen fl ° W s * tua t ion as a 

^' 3 — e Modeling Principles and Asstm^tjons 

kl the derived mean momentum and 

kine ,c energy conservation equations, i t can be 

“ s be , th3t SOr " "O d *“'« assumptions 

a7d 7h ' on the physical interpretation 

arjd the nature of each term. Previous 

xperimentai and theoretical findings can help in 

^ "f, C ° UeCti r el y" ‘-ms with minimum 

e m P lrlc al constants. Secondly, carrying 
out an order of magnitude ar^lysis for diffe^M 

equa t ions° nS apPW ln the S°vern^ 

th^mt S T termS ray ** neglected. Thirdly 
the micromechanics which control the ability of 

S the° W f; a f ri3ble rr t0 " lth other 

correlnrt C °7 af(>ctIn « the magnitude of these 
correlations should be considered. With the 

^ of S si r 17i rkS ln K mlnd - ^ “» «•« the 

saj<e oi simplicity that: 
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1 ) 


The correlations between the Urge eddies 
from the production range and the snail 
eddies from the transfer range (mixed 

correlations e.g. uu 77 ) can be neglected as 
they originate differently and they are 
related to different ranges of the power 
spectrum. Similar assumptions are accepted 
in the classical single fluid turbulence 
theory. 


where 


Yp = C mP p L k LP 

^LT = C mT P L k LT 

P x T 

•Yt = ‘Yp + Yt 


1/2 1 
V 

(15) 

1/2 , , 
l.j. . and 

(16) 

= C u P L k L 1/2 >L 

(17) 


2) The void fraction fluctuations occur mainly 
at low frequencies i.e. 


a = a + a. + a 


3 ) 


( 10 ) 

( 11 ) 


with a’ >> cr M 

This Is a simplifying assumption which is 
acceptable for such complicated problems. If 
the particles are of small diameter their 
concentration is relatively uniform 

distributed in the Taylor length scale. High 
void fractions are mostly associated with 
large size and high density particulate 
flows. These large size particles are mainly 
fluctuating at low frequencies due to Its 
high inertia, hence they in turn correlate 
weakly at high frequencies. 

The correlations of higher order than three, 
for Instance au'dp'/dx 


i ’ 


• ♦ 

u .u 

i j a* 


etc., are neglected. These are at least an 
order rmgnitude smaller than those of the 
third order (see Hanjalic’ and Launder 
(1972)). 

4) Pressure diffusion contribution to the total 
turbulent diffusion in the kinetic energy 
equation will be neglected because of its 
relatively small magnitude (Hanjalic and 
Launder ( 1972) ) . 

5) The Bousslnesq gradient type approximation is 
adopted for modeling of different fluxes and 
triple correlations, with assumptions similar 
to Elghobashi and Abou-Arab (19S3) and Roco 
and Mahadevan (1986) 

6) The following constitutive relations are 
employed for the shear stress of carrier 
fluid: 

- = Yp ( 55T + ' 3 k LP*ij" 3 YphjVn 

j (12) 


- pu lUj 


chi du r, 2 — 

^dT + ' 3 k LT fi ij ■ 3 YT 6 ij U n.n 

J (13) 


while the total shear stress -p u.u^ can be given 
by 

(14) 


Similar relations can be written for the viscosity 
of the dispersed phase (see Roco and 

Balakrishnan (1985)). However, in the present 
work we choose to define the eddy viscosity of 
solids as follows: 


u St = u Lt 7 ° 


where 
and 


° = a oS 7 °S 


°aS = 




( 18 ) 

(19) 

( 20 ) 
( 21 ) 


Appropriate expressions for are cited in many 

articles such as Peskin (1971). Picart et al. 

(1986). and Hetsroni (1982). cr^ is a Schmidt 

number and its value is about 1.5 (Abou-Ellail and 
Abou-Arab (1985)). 

7) In the present approach for dilute 

particulate flows the turbulence kinetic 
energy equation Is written only for the 
carrier flow. The solid phase turbulence 
kinetic energy and turbulence correlations 
are evaluated from the available solution of 
the linearized equation of motion of a solid 
particle after its transformation from the 
real time to the frequency domain. 

8) Terms which are of similar nature i.e. 

convection, diffusion, dissipation, etc. can 
be modeled collectively. The length-, 

velocity— and time-scale which are 
appropriate for the description of their rate 
of change must be identified from the 
physical interpretation of these terms. 

9) The response function which shows the ability 
of solid particles to follow the eddies is 
obtained from the equation of motion of 
particles for different local dimensionless 
parameter, dg/1^ and p^/p^. 


10) To establish the degree of generality of the 
proposed model validation tests were carried 
out for air-laden and water-laden jet. and 
air-solid pipe flow. 


5. CLOSURE FOR THE MEAN FLOW EQUATIONS 


With the modeling assumptions given in the 
previous section, the steady state mean flow 
momentum equations for any flow component, reads 


P K <3x, ^ a K U Ki U Kj^ + dx P K (VYiYj * a K U Ki U K j 


Uv.) 
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+ ^ P K [(Q X u K, u Kj + a K u Ki u Kj^ + u Ki (a K u Kj + “JMcj) 
+ U Kj <VVi + %)] + 

= P K ( °K b Ki + “K^j + a K b Kl )~ ^“kPr 

* (F k ^ + Pk ^ + ^ S) 


The diffusion fluxes namely tyi^ . a^uj^ etc. 

are modeled using Bousstnesq approxiration as: 


— — = \? 

Ki p &x . 


while 


°aK 


— — ^ITT 

°K U Ki = ~ fa 


°aK '" i 


(26) 


(27) 


9x (a K T faji ) + (I Ki } -K 


( 22 ) 


In Eq. (22) there are 16 correlations, half 
of them in the production "large eddy" range of 
the spectrum The terms are modeled following the 
criteria. 1 ) Physically correct behavior, ii) 
Minimum number of empirical functions and 
constants. and ill) Comparisons against 
experimental data over a wide range of conditions 
are required to check the validity of the model. 

The turbulent stresses caused by the large 
and small size energetic eddies. and 

~ P L U i u j- are defi ned by Eqs. (12) and (13). The 

correlation between uj’ and uj 1 is weak when i / j. 

This finding will be explained in the next 
section. 

, Tbe colllsional effect correlations are 
modeled after Launder (1976): 

‘ C .5 <=kS 7 3^ =K?X 


U k'j u Kl dxj 


(23) 


where is a constant of a value of about 0.1 


Simi lar ly 


,Sa 


a K U Ki u Kj ~ ' C *S ( t KT ) (u Ki u Kl 


Similar expression can be written for ^ and 

V'KJ etc - However. according to our 

assumptions, especially those concerning the void 
fraction fluctuation at high frequency a", the 
fluxes in Eqs. (26) and (27) can be modeled 
collectively as : 


f3< 

^ Kl °oK ^i 


(28) 


,, ^ sal “ tl ° n of the two transport equations 

* ° r k p “"d kj) would require also a description 

for the length scales. This point will be 

explained in more details in the next section. 

The fourth group of terms with the void 
fract on-shcar stress turbulent correlations 
" v the lajninar viscosity as a multiplier 

and will be neglected due to its smaller order of 
magTii tude. 

The pressure effect contribution to the mean 
flow equation consists of two groups. Each 
contains three terms. The first of these terms Is 

rh TTT P ressurc -v°id fraction. The second and 
the third terms are the pressure-void fraction 

and r !^^ Hons which can be modeled after Elghobashi 
and Abou-Arab (1983) by 


" (a K p K * = *1 + * 2 

where 


(29) 


= - C 

and 


1/2 


* n l, 1 ' t t 

*3 P K \ ‘ u Km a K 


+ U Kj u Kl 


U Ki°K> 


( 24 ) 


*2 = " Vk k 


1/2 


t t 


where C^ 5 must be optimized by comparison 

experimental data The fluxes in (23) and 
can also be collectively modeled as: 


wi th 
( 2 - 1 ) 



The values of the constants C and C are about 
unity. * +* 


me second correlation in the second group of 
terms can be also modeled following Launder (1976) 
Ine final form is J 


^S +p i^=P K <^K CV'W 
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t t 
u.. .u 


^ , Ki Kl 2 * ,. M t „ 

+ C *2 ( " 3 6 il ) K1°K ] 


♦ ^ 8 « • ' °' 2 “M • > (3 ° } 


The values of the constants and C^ 2 are 4.3 and 

-3.2, respectively. The modeling of these 
correlations suffers from the embodied assumptions 
concerning the velocity and length scale 
description. It would require a large number of 
transport equations to model accurately each of the 
above correlations. 


The interaction term ( I^i )-K for * = ^ 

(liquid) and -K = S (solid) is modeled for dilute 
suspensions with particle Reynolds numbers less 


suspens 
than uni ty - 


,2, r/ - 


(I^ = - (MV^) - u Si^ a s 


,\p + t f!k 

+ ( P T * dx 
°aS °aS 1 


aa aS °°aS ‘ 


where the first term is the drag Interaction for 

particles in the Stokes range. The second and 

third terms are the turbulent fluxes due to the 

relative motion between the particles and fluid. 
The gradient transport model with the exchange 
coefficients \> L? and corresponding to the 

production and transfer ranges is adopted for 

these fluxes ( a£ u^ , u l,- a s u Si' 

a S U Si } - 

If only single velocity scale is chosen for 
the whole energy spectrum, k^, there will be only 

change coefficient instead of 


one 

momentum 

U LP 

and v LJ 

V Lt 

- VL? + 

P 


oo „ 
aS 


LT 


J oS 


(32) 


6. CLOSURE FOR THE KINETIC ENERGY EQUATION 


The first group of terms in the k-equntion 
(Group ttl) is the convec t ion of the total specific 
kinetic energy, where 


ak L 

dZ 


ak 


LP 


j 


ax . 

j 


ak 

dx 


LI 

j 


(33) 


Pi f fusion transport of k is composed of two 
main parts. The first part (Group *42) is the 

velocity diffusion and it contains the 3rd order 
velocity correlations, while the second part (Croup 
t*3) is the pressure diffusion with the 
pressure-velocity correlations. The modeling of 






_ - U LP ^LP 

P L “L u Lj u Li u Li = ' P L°L ~P~ dx 

°k J 

and 

_ - U LT 

P L °L U LJ U L1 U L1 = ' P L“L ~r 

°k J 


(34) 


(35) 


where a?* and oT are the turbulent Prandtl /Schmidt 
k k 

numbers for the kinetic energy in the production 
and transfer ranges. To reduce the number of 
empirical constants and the number of governing 
equations the above two correlations are modeled 
collectively as follows 


- P Lt 


ak, 


~ p L U Lj U Li U Li = - p L°L 


a. dx . 

k J 


(36) 


where o. is of order of unity, 

k 

The presure diffusion term is negligible 
relative to the velocity diffusion (Hanjalic' and 
Launder (1972)). 

Mixed and higher order correlations (Croups **3 
and #5) can be neglected according to the modeling 
assumptions stated and discussed in section i.3 of 
the present paper. 

The production terms are divided into two 
groups. The first group (Group t*6) is common for 
single and multiphase incompressible flows. 


du 


Li 


au 


Li 


P L°L u Li u Lj • ET~ * P L°L u Li u Lj * ftT 

3 

be modeled collectively as follows: 


and can 


In the present work, the turbulence kinetic 
energy for the liquid phase "k^” (turbulence 

velocity scale) is obtained from an exact transport 
equation (Eq. AS in Appendix A), and the length 
scale "l** is described algebraically. The kinetic 
energy equation AS contains a large number of 
turbulence correlations. In order to obtain an 
engineering turbulence model, it is sufficient to 
consider the principles and the assumptions given 
in Section 4.3. By engineering turbulence model, 
it is meant, a physically correct model with 
minimum number of empirical coefficients. 


du. 


du 


P L°L u Ll u Lj ‘ 0 L m Li ( 5x 


Li x 2 


(37) 


The physics of turbulence and the 
consideration of the spectral energy transfer 
assume that the production is only due to the 
interaction between the mean flow and the large 

eddy. Since the u^.u”~ correlation is for medium 

size eddies which have almost no direct interaction 
with the mean flow. it results that the 
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contribution of these small eddies to the 
turbulence production via the mean field is smaller 
than that of the large eddies. This 

correlation is weak If i ? 


U L1 U LJ 


means also that 
J. !•«- it is only of significant value if the 

o 

turbulent normal stress components (uj^ . 1 = 1, 2, 

3) are considered. According to Kanjalic' et al. 
(1979) multiple scale model the turbulent viscosity 
is defined as follows^ 


= S k L ( k L/ £ Lp) 

(38) 

1 = k LP + k LT 

(39a) 

= 0.09 

(39b) 


and C 


This equation can be rewritten as 


du 


like ^ and o£u£ ( ~ — are modeled In 

a similar manner to that of the above terms i.e. 

- (°lT 7 + aj'u" ) — ii = ^ ^ 

J Ej 3x o dx . dx . 


(45) 


^ ~ ra Production terms (Group #7) can be 
written In the following form 1 can be 


Extra production = -tl ( fllk + “ ^Li 

°aS dx , LJ aT~ ) 


°aS U fttj } 


(46) 


^t^L 


C (^ , S r > 

p £ lp £ lp 


= C pP k LP *LP + C p1 J< LT 5l LT = ^*Yp + ^T )/p L ^ 

where C^p and are two additional constants and 

*LP *LT arc a ^ so two additional length scales 

for the large and medium size eddies. The length 
scale can be related with the following relations 
(from Eqs. (3 8) and (40)) 


'lT " C 1 J LP * k LT /k Lp) 


(41) 


Since the ratio I s of the order of unity 

*LP > *LT* the constajU c j should be smaller 

than unity. Thus If the multiple time scale model 
is not recommended (due to its large number of 
additional constants) an alternative approach is to 
consider a multiple velocity scale model. In this 
model only two differential equations for k Lp and 

k LT solved. The length scales can be 

obtained by using the previous relation Eq . (41) 

and any expression for the length scale of the 
large eddies l^p, for example that used by Roco and 

Shook (1983) for cylindrical pipes. 

The model ing of the additional production terms 
(Group #7 ) is achieved as follows : 


*P, 


LP 


- <3a L 

°L U L1 dx t ~ T &T ET (.V) 

°aS 1 1 


. ~ — - — U LT ^ 

“i-^l dx ~ T dx dx ( 43 ) 

1 CT aS 1 1 

and both collectively as 

. . ®Pi u i . 3p, 3a. 

<«) 

1 aS i i 


_ chj.' . du 

The terms p, cl ll — ~~ Li „ — 

Wl ax ax W L §7- 

J j 


du" du". 
Li Li 


in 


Croup #S represents the dissipation rate from large 
^dtes e Lp and transfer ^Tnce viscous 

dissipation Is mainly confined to small scale 
~?!ini an M t0 S,n T ify the >®thematical form of the 
al ^1979)"* h! 3 * turbulence m °3 e 1 s . Han ja lie' et 
eoulllhM 79) haVC assumed th a' there is an 
dlsi Sl ? c,run ’ 'Clergy transfer between the 

dissipation and transfer region i.e.. t = t 

where £j _ i s the dissipation rate in 'the* ' s i ng 
scale scheme. Thus 


and 


p lV 


chi,’ ^ 

du 1 

Li 

Li 

dx. 

dx A 

J 

J 

au" 

du" 

Li 

a. . 

Li 


a5 


Li_j2 

J 


l ax dx ~ Wit 


(4S) 


and 


where 


£ LP = Sp k Lp 5/1 LP 

fc ET - Sr * 5 - lt 

JrodTuon" anJ ) 

(Han jai i c ' et al f 197911 n“ . Considered 

’ IP : q e r™ ~ 

veloc^v P ° rrelations between the fluctuating 

forces (intera°M nt ^ the fluctuatln £ friction 
due to V f Tr. r trrm - ln Crou P s ** and #10) are 

dHute flow's ^ T d f / Uld ‘ S ° lid 1" 

to t frlctl on interaction terms due 

CrouT ^ Coll,slon (fluid-fluid interact ton 
Croup HS ) are given above by Eq (47) and (4S) 

The form of the correlation between the fluctuating 
ag force and the velocity fluctuations depends on 
the expression adopted for the drag force The 
viscous drag correlation (VDC) in Croup #10 J^ r 

rea5s S ^ in dilute suspensions 
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VDC 


■ - ^ l“u - =».>•*. 


with the following dimensionless parameters 


s = P S /P L 


18 ^L 


(“S' A 2 + A 3 } 


s 

where** 


A 1 = °S U i + “S U 1 


a 2 = , (u. 


LI ' Ll “Si 


u c.) + u Li< u Li' u Si ) 


(49a) 


(49b) 

(49c) 


*S = 


l>^ »c I xRe 


“T 


A 3 = a S u Li (u Ll" U Si ) + a S U Ll (U i'l~ U Sl^ (49d) 

The first correlation group Aj (Eq. 49a) can be 

approximated using the gradient type assu^tion. 
The second correlation in this expression (49b) is 
that due to the relative slip fluctuating motion 

u Li( U Li - U s7> “ d U Li (u Li * U Si^’ 17,656 ^ ^ 
modeled using similar approach to that of 
Elghobashi and Abou-Arab (1983) but with some 
modifications which allow for different 
particle-eddy interaction according to their 
relative size. These modifications are based on 
the spectral analysis carried out by Dingguo (198f) 
for the response of the particles to the turbulent 
fluctuations of the carrier fluid. For very large 
eddies k << 


where 

= is the wavenumber of the most energetic 
eddies. 

Re = Reynolds number based on l^, and 

N = Stokes number. 

S 

For small eddies with wave numbers x >> x g . 
the particle response can be described by 


k S,3 P L 


p; 


(53) 


For the intermediate size eddies * - *5 one 

can use (SO) wi th 

1 


a = - and 
s 

0 = tan 1 (s) 


(54a) 

(54b) 


as 


If the fluctuating slip velocity w j is defined 

(bo) 


= U S> - U Li 


u g = (v s ) K = (v l ) k a • exp[-i(xu L t-p)] 

(50) 

then the 

ratio 

of 

the mean 

.2 

square w 

a ■ 2 

and u L . 

where k is the wave number; K g 

is the 

becomes 






Basset-Boussinesq-Ossen wave number defined as 
l/d g : ( v g) and (v^) k are solid and liquid 

velocity components with the wave number x: a is 


= 



u Li> / U Li 

(56) 

“2/~2 
w i /u Li 

- 

2 h^si + 

the amplitude ratio of oscillations 

,2 2.0.5 

a = [{1 ♦ q 2 ) + dil 

(51a) 

with 







1 

— 

~~ 9 . 

u Si 2/u Li 2 ' 

"7 2 ~ 2. 

w i /u Li> 


u Li U Si 

= 2 

U u V * 


and P is the phase angle of oscillation 

P = tg _1 [ci 2 /( 1 + dj)] 

The expressions of q^ and q£ are 
9N S w 1-s , 

(1 + H s+ o 5 

-l2( s+0.5) s 0 5 


(51b) 


= k U Li ^ + r SL " r RL^ 


(57) 


61 i2 (^s + |) 2+ (« + 


9N 


(52a) 


The values of r g] _ and can be obtained by using 

the preceding solution (Eq. 50) 

r gL =[| S a 2 E L (K)dx + f (-f) 6 s'\( K ) d ’ t ) / J^ E L( K)dK 


S ) 2 


(s+0.5) 


-f2( s+0 . 5) 




or 

r 


g( l ~ 5 )^* ( 2 N? + 5 ^) 


"2 T 2 
( s+0 . 5) 

. fi 


9K 


8 i_ (2N 2 + V2 __2S_ } 2 

0.51 2 S J2 T2( s+0.5) 


(52b) 


SL ” w(h+ 1 ) 
where 


2h tan '(Kg/x-j.) 


(58) 

(59) 


(s+0.5) 


h = 


■c-jU^ (s+0.5) 


18 (xg/x^/Re. A = 


lSu, 


(s 4 0.5)d^ 
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Figure 2. The effect of eddy-particle size ratio 
'VV on the Particle response to the 

n „ Q 

eddy fluctuations (u^/u^) 


The above analysis applies equally well to the 
large eddies as to small eddies. The energy 
spectrum function E^k) for a two-phase flow is 

given by Al-Taweel and Landau (1977). However, 
since its form is not essential (Dingguo (1967)) it 
is sufficient to adopt any simple form as that 
given above by Eq. (61). It is clear from the 
above analysis that the particle response to the 
carrier fluid fluctuations is a function of the 
density ratio Pg/p^. size of interacting eddy 

relative to particle size k^k and Reynolds number 

based on the size of the most energetic eddies of 
the flow. An analogous expression to that given by 
Lq . (62) is that based on the Chao's solution (see 
Chao (1964)). This solution can be considered as a 
substitution of Dingguo's solution only for k << 

i e. for fine particles. 

The last term to be modeled in the VCD group. 
^ 3 . is separated into four correlations' 


^ 

h = - -J- [f“s u U u L> + a s\iV;> 
s 

T 1 T2 

' (a S u Li u Si + a S u Li u Si^ 


(63) 


and 

r RL = (a 2 l) E l («) d * 

+ r / r £(*) dK (60) 

K S J ° 

where = (f^ + f^) and E^(k) is the liquid 

energy spectrum function for which any arbitrary 
form can be adopted, for instance 




.2 

2 u Li 


*T 1 




(61) 


Figure 2 illustrates the effect of 
eddy-particle size ratio expressed as on the 

2 2 

ratio f SL~ U Si ^ U L\ w ^ich indicates the particle 

response to the eddying motion. It can 

be noticed that high values of (i.e. small 

particle or large eddy) the particles follows quite 
well the eddy motion. 

Substituting the expressions for and f 

(Eqs. (58) or (59), and Eq. (60)) Into Eq . ( 57 ), 

the correlation can be obtained 

A 2 = l u U (1 - f SL + r RL ) («) 


T3 


T4 


where the triple correlations T3 and T 4 are modeled 
in a similar manner to that used for the 
calculation of Tl and T2 by using Eqs. (23) and 
(24). thus 


Tl 


^LiV 


°S U L 1 U L 1 = ^5 (“LP^LP) LI P- 




72 = “sWi = -Ss -gf - ) T - 


T3 = -C 




Vs (k LP /t LP )(u Ll u Si — 


♦ u’.u 


•*Ll a S . 
'Si u Ll —a^~ h 


T4= -c 


>5 {k LT /t LT )(u Ll u Sl 




dx. 


■ + u" . u 


^Ll a S 


Si LI d^~ >1 
(64) 


These correlations can also be collectively modeled 
using single velocity and length-scale, and total 
L and £ L‘ The first two and the last two terms 
yield, respectively: 


2 k. 


T . . 72 . c>5 _t 


L 

2 k, 


>5 


, t t 
(u Li u Ll 

Li ! 

*1 


o' 

(u Li U Si 

ax, 


(65a) 
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i.o 


u Si U Li 


- t t 
dx . 


(65b) 


The terms in Group #9 of Erj. A5 are of 
diffusive and dissipative nature. The diffusion 
terms as they appeared in Group #9 are multiplied 
by the molecular viscosity and therefore will be 
neglected due to their relatively small magnitude^ 
Other higher order correlations and nnxed 

correlations in this group are also neglected 
according to the modeling principles state 

prevloulsy in section 4.3 and as they are also 
multiplied by the molecular viscosity. 

By substituting all previously modeled terms 
into the exact form of the turbulence kinetic 
energy equation, and rearranging these terms, one 
obtains the simplified modeled form given in 
Appendix B (Eq. B2) . 

Since the present model is based on an exact 
equation, namely the turbulence kinetic energy 
equation, and the modeled form of this equation has 
no adjusting coefficients it is expected that the 
model will generally produce good results and have 
less limitations compared to other models. The 
only modeling assumption is that of the Boussinesq 
gradient type, which generally is accepted. The 
correlations that requires questionable 

semi-empirical modeling assumptions and 

introduction of empirical constant are (i) fewer in 

number (for a L U S i U L j ^ ^° r 

terms having an order of magnitude smaller (by 

ratio a' /a) compared to other main terms in the 
k-equation. The only significant new correlation 
used in the present closure is that due to the 
relative motion between the phases. This is 
modeled with less restrictions and taking into 
consideration the effect of the particle 
diameter-eddy size ratio on the particle response 
to the eddying motion. The limitation of the 
present one-equation k model closure is the 
algebraic formulation for the length scale. Since 
there are many factors affecting this length scale 
and since it Is even difficult in many practical 
applications to give a unique and accurate 
description of the length scale, the use of a 
transport equation for the length scale in the 
two-equation model of Elghobashi and Abou-Arab 
(1983) is expected to give better results with fine 
particles (d g < q) . However, it should be noticed 

that the former model and any other similar models 
contain some empirical constants specific for 
various flow conditions, and they require the 
solution for an additional transport equation. It 
can be expected that the present model combined 
with an appropriate length scale equation e.g. 
dissipation rate equation will simulate better most 
of the important . features of multiphase turbulent 
flows, particularly the fluid particle interaction. 
In that case the number of closure transport 
equations will increase to three (if two velocity 
scale, kp and k^.. and one length scale transport 

equations) or four transport equations (if two 
velocity scales and two length scale. Ip and lj. 

are adopted). 



Figure 3. Computed and measured [20] mean velocity 
distribution of air in air-solid pipe 

flow ( — •• single-phase. 2 Eqs. k-e 

model [12], — 1 Eq. k-model [27]. -xx- 
1 Eq. k-model [12]. 1 Eq. i> t -®odel 



U L 



Figure 4. Predicted radial distribution of the 

turbulence intensity of air in air-solid 
pipe flow using different turbulence 
models (keynote as in Fig. 3). 


7. SAXPLE OF RESULTS 

Figures 3 to 6 compare the present predictions 
with LDA-measu remen ts for single and two-phase 
turbulent pipe flow (see Haeda et al (I960)) and 
turbulent round water Jet laden with uniform-size 
solid particles (see Par thasarathy and Faeth 
(1987)). Both flows are oriented vertically 

downward. These flows are axi synrotr leal . The 

concentration profiles are given as Input data 
based on experimental results. 

In these flow situations the average eddy size 

£ was varying from about one d c to few hundreds 
e 
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d g- The corresponding representative eddy size in 

the transfer range was only a fraction of particle 
diameter in the pipe flow case. In the jet flow 

case the mean size of large eddies and the 
Kolmogrov length scale were also varied in a wide 
range. These scales are field variables, and they 
depend upon the flow configuration, location in the 
flow domain and particle dimensions. 

Two-phase flow solutions were obtained by 
solving the flow governing equations in their 
modeled form which are described in the previous 
section and given in Appendix B. The numerical 
procedure used for these predictions is based on a 
developed version of the Genmix-Code of Spalding 
(1977). However, since the main objective of this 
paper is to give a complete description of a 
developed turbulence closure for multiphase flows 
and due to the space limitation the details of this 
numerical approach will not be given here. The CPU 
Time for the two considered flow cases was about 4 
minutes on a VAX 760 Mini-Computer and 2 minutes on 


I- Gas-Solid Vertical Pipe Flow: 

Figure 3 shows a comparison between the 
experimental data and the present predictions using 
five different models of turbulence namely 1. 
One-equation k-model of Roco and Kahadevan (1986). 
2. One-equation u ^ -mode 1 of Roco and Ba lakr i shnan 

(1985), 3. The k-equation as given in the 
two-equation k-t model of Elghobashi and Abou-Arab 
(1983). 4. The two-equation model of Elghobashi and 
Abou-Arab (19S3). and 5. The 
present one— equation k-model. The figure displays 
the mean axial velocity distribution In the fully 
developed zone of the pipe flow for single and 
two-phase cases. The differences between the 
predictions of all one-equation turbulence models 
and experiments Is mainly caused by the general 
algebraic expression adopted for the turbulence 
length scale which was not optimized or adjusted. 
In the present computation the concentration 
profiles are assumed based on previous experimental 
data. The inlet concentration distribution Is 
taken to be similar to that given by the best curve 
fit after the experimental data of Soo (1967). 
Figure 4 comperes the calculated turbulence 


intensity defined as J\ 2 (=§ \) / u L with Its 

measured values. The near wall treatment is based 
on a modified form for the law of wall (see 
Abou-Ellail and Abou-Arab (1984). Lee and Chung 
(1987)) and the particle slip condition at the pipe 
wall. 


The present model as compared with the one 
equation models of Refs. [12] and [27] predicts 
slightly higher values for the mean flow 
quantities. However, it must be mentioned that the 
last k-models contain empirical constants in the 
dissipation term of the k equation. Concerning the 
fluctuating flow quantities the present model gives 
slightly better results for the turbulence 
intensity in the near wall region than other 
one-equation models. The expression for the 
turbulence length scale was not optimized. It Is 
also expected that the current model will give 
better predictions for coarse (d c > tj) and heavy 


U L / U 



r / z 


Figure 5. Computed and measured [23] mean velocity 

Ji M f 1 r ftrxw 



r / z 


Figure 6. Predicted radial distribution of the 

turlxilence intensity of water in water- 
solid Jet flow. 


suspension flows for which no set of comprehensive 
2D data is available for comparison. 

Case II : Turbulent Round Water Jet Laden with 

Uniform-size Solid Particles 

Comparison between experimental data and 

numerical predictions of different turbulence 

closures are given in Figures 5 and 6. 

The two-phase flow measurements on velocity, 
concentration and turbulence correlations used for 
comparisons are taken from Par thasarathy & Faeth 
(1987). Different axial locations within the round 
Jet, between eight and fourty Jet diameters from 

injection nozzle are 
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r / 2 

Figure 7. Radial distribution of solid phase 

concentration in water-solid Jet flow. 


considered. The radial distribution of solid 
concentration given in Figure 7 are at eight 
diameters from the nozzle. 



0.2*10'* 10'* I0' 3 I0‘ 2 I0’ 1 I 10 d s // (5*) 

(b) 


The comparison shows that the two-equation 
two— phase k-t model of Elghobashi and Abou-Arab 
(1983) describes both flows better than the 
one-equation mass/time averaged turbulence model. 
However, it is Important to mention here that the 
presently developed closure uses only one transport 
equation and without adjusting any empirical 
coefficient. At the same time, the present closure 
is in its early stages and more refinements and 
validation tests are required. especially for 
coarse particles two-phase flows for which one 
would expect that the present k-formulat ion will 
provide improved predictions. The difference 
between the predictions of the mean and fluctuating 
flow velocity components as obtained by the present 
k- f ormu 1 a t 1 on and that of Ref. [12] depends on the 
particle size relative to its surrounding eddies. 
Figure 8 illustrates this difference at two 

different loading ratios for the pipe flow at a 
radial distance (R - r)/R equal 0.1. 

6. CONCLUDING REMARKS 

The turbulence closure 
for dilute suspension flow 

turbulence kinetic energy 

features of this model are 

i) Two velocity scales are adopted in computation 
for large and medium size eddies, 

corresponding to the turbulence production and 
transfer range. respectively They are 

expressed into the governing equations by a 
specific local mass/time averaging On this 
basis the spatial and temporal transfer rates 
of the thermodynami c quantities and the 
particle-eddy interaction are better 
es t imated . 


Figure 8. Difference* between the present 
k-formilat ion and that of Ref . 

[12] for different mass loading ratios 
L.R.: (a) fluctuating and (b) mean 

axial velocity. 


ii) Spectral analysis of the interaction mechanism 
between particles and most energetic eddies 
provide analytical correlations for closure 
The particle response and the modulation of 
turbulent eddying motion is given as a 
function of the p>ar t ic le-fluid density and 
size ratios. 

lii) To keep the number of transport equations of 
the turbulence closure ar.d the number of 

empirical constants as minimum as possible, 
the length scales lp and 1-j- are described 

using algebraic expressions. Relations 

between these scales (Eqs. 40 and 47) are 
suggested . 

iv) The model does not introduce additional 

empirical constants to the closure of the 
velocity scale equation. 


presented in this paper 
is based on the fluid 
equation. The main 
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where 

K = phase (or generally a flow component) 

T = mass/time average of f over K 
K 

i.j = 1,2,3 (Cartesian coordiantes) 

H,. = body force in the i-th direction 
K i 

(T~~) ^ projection in the i-th direction of the 
K i — K 
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APPENDIX - A 


interaction vector 


The interaction term as it stands for 
sol id/1 iquid drag is given by 

Ss l u Lv‘ U Sv^ U Ly' U Sv ) 

(I sA - "u’s ■ o -"Vlt; ' ( ‘, ,v:t— 


Mass/Time Averaged conservation Equation s 

Mass averaging the conservation equations of mass 
and momentum over- a flow component (K) one obtains 
a new system of equations for mean velocity, 
concentration and kinetic energy of turbulence. 

The point Instantaneous conservation equation 
can be written for any flow component (K) or for 
the entire mixture in the folowing general form 

(p4<) ♦ v(p'J'u) ♦ v*J - S * 0 (Al) 


(A3) 

This drag term takes a simple form for Reynolds 
numbers less than unity 

(I Si } L = ' {I Ll } S = “ s (l 8 H 1 / d s )(u Li' U Si ) (M) 

The transverse effects caused by the presence of 
other solid particles. Saffman force and Ho and 
Leal inertial force are neglected. Their 

Importance is small in dilute suspension turbulent 
flows with fine particles. 


where : p = densl ty 

u = velocity vector 

= transported quantity 
J = flux vector for ^ 

S = source term 

Let assume * * By splitting each flow 

property into mean and turbulent fluctuating 
component (u^ ♦ u^) and mass/time or double time 

averaging the equation (Al). one obtains the 
following momentum equations in the i-th direction 
for an incompressible phase (K) without mass 
exchange with other flow components (see Roco and 
Shook ( 1985) ) : 

P K It VKi> + p K (“K“Ki' r Kj> 

Time Rate Mean Flow Convection 


t -- _ t da^ 

=p K°K^Ki ■ dT x (wW * h ftT + Pk ftT - 

Body Force Pressure Effect 


* dfj ~ p K a K* U Ki U Kj 

Frictional Inertial 

Effect Effect 


Equation (A2) contains terms due to the 
unsteady flow, mean flow convection, diffusion, 
pressure. body force, as well as frictional, 
inertial and collisional effects. The mean form of 
the turbulence kinetic energy governing equation is 
obtained by subtracting from the steady state 
instantaneous momentum equation for a component K 
L the corresponding mean flow equations, and t en 
multiplying the resulting difference equation by 
(u^. + u**^.). By averaging one obtain the kinet \ c 

energy equation for a flow component K - L. This 
equation reads 


- - ,^Ki /2 . 

p K a K u K J ~dx^ + 

Group #1 (Convection) 


♦ 



Group #2 (Velocity 

Di f fusion) 


+ U Kj U Ki U Ki 


+ u kj u K> u ki 


u Kj u Ki U Ki 


u,, ,U t , . u. . 

hj ki hi 


+ u Kj u Ki U Ki + u Kj U Ki u Ki + U Kj U Ki U Ki )] 

+ (WKi U Ki u Kj + a K u Ki U Kl U KJ ) 

Croup #3 (Higher Order Correlations) 
+ other 4th order and minor terms) ]= 
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°K dXj (u Ki P K + \l ? K + u Ki P K + u Kl P K ) 
Group tH (Pressure Diffusion) 


_ Du 
°K ( P K dx 


w . Du” . Du” . Du ' . 

Kl + p;;^ 1 + p--# 1 + p;> Kl 


K dx J 


K Dx, 


K Dx. 


i " i " i i 

Group ttS (Extra Production and Transfer) 


. ^ 




^“k d* i * dx t + Sii'Si alj 


^p k <?p k 

+ V^i dx~ + I*} 


9u 

(p K a K u Ki u Kj ~ * p K a K u Ki u Kj aT~> 
J J 

Group #6 (Production) 

r- ^ ^Ki — T — 

' l u Ki “k u kj &r~ + “kjVki-ST 

J J 

Croup #7 (Extra Production) 

_ __ Du^ . 

* “Ki-lMcj “*rr + “Kj'SMciTST 

3P K 1 

" (a K u Ki + a K u Ki ) dx - + lninor terms I 


(t - - ^ia Da 

Wk ST ~ET + Wk — ~ 
J J J J 

Group #8 (Dissipation) 


Di Di 


Du” Du” 
X: Ki 


+ vvk sr nr * vSD -§t -*r + mi « d 

J J j j correla- 


Group #9 (Extra Dissipation & Diffusion) 


t ions 


^iDd-K + SaDD-K 

Group #10 (Extra Dissipation) 


(A5) 


APPENDIX B 


i ne Modeled form of 

Enert 


_the Turbulence Kinetic 
mat ion 


The steady-state turbulence kinetic energy equation 
for the liquid phase (K = L) is: 

- - D D ‘‘Lt D 2 

Wlj KJ-RJT-KJ* Vl, (s ^-) 2 

Convection Diffusion Production 


♦ “ti fk { !!k + G ^Ll U Lt D - 

°aS °*i lax , U (U U 

Extra Production 



Dx , 
J 


) 


Wt 

Di ssipat ion 

18^ 

~2 


"Lt D 

“S ~' aS 

Extra Dissipation 


Di _u Sib 3x. 


18a, 

■ VD^SL^RL* - ~Y- (T1 ♦ 12 - T3 - 

d s 

where the expressions for J\, . r . T1 . 

oL KL 

T'i are given in the text. 


T4) (Bl) 

T2 . T3 and 
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